In many of the predictive and analytical problems of theoretical and applied engineering, it is necessary to solve one or more of the equations of physics that describe a physical field. The analyst will target the solution of one or more field variables, for example, mechanical deformation, acoustic pressure, electrostatic or electromagnetic potential, or electromagnetic field intensity, among others.
The equations that govern these and other fields are differential equations, integral equations, or equations that combine aspects of both, which can be solved exactly, if at all, only for specific, simple cases. Boundary conditions (and initial conditions, if time evolution is part of the inquiry) complete the definition of the problem. In practical situations, these boundary conditions often involve specifying the field variable, its derivatives, functions of either or both of these, or a combination of some or all of the preceding, on complicated two- or three-dimensional surfaces. Due to the complexity of both the governing physical equations and the boundary conditions, most practical problems are too complex to solve by hand calculations.
Numerical modeling techniques are performed with the assistance of a digital computer to provide approximate solutions to these equations. Many different numerical modeling techniques are available. One such technique is the finite element method (FEM).
As we describe more fully hereafter, the traditional practice of the FEM calls for a spatial region of interest to be subdivided into a plurality of perfectly juxtaposed cells (i.e. cells that are juxtaposed without overlaps or interstitial voids), referred to as elements. Discrete nodes, where solution data is tracked, are defined on the inter-element boundaries and, sometimes, also within the element interiors. We refer to the network of elements and nodes as a mesh.
The FEM does not seek an exact solution to the physical field equation. The FEM assumes that within each element, the field can be described adequately by a summation of simple functions, such as polynomials, that are chosen, inter alia, for both convenient analytical properties and to provide an adequate approximation. The coefficients of the respective polynomials, corresponding to governing solution quantities at the nodes, in such a summation are referred to as degrees of freedom (DOF). The polynomials are frequently normalized in such a way that the DOF are equal to, or otherwise relate to, the values of the field variable at the nodes. Between the nodes, the polynomials interpolate the values of the field variable in a continuous fashion.
The polynomials themselves are known a priori. Therefore, the field values within a given element are determined by specifying the values of the DOF for that element. In FEM, the DOF values are not determined directly using the field equation; instead, a mathematical condition is derived from the field equation which imposes that a certain measure of error (i.e., between the approximate and exact solutions) must be small. Applying this condition to the field equations results in a system of equations in which the unknowns are the DOF. These equations are well suited for solution by a digital computer, since they can be solved via s a large number of repetitive, mechanical manipulations of stored quantities.
The system of equations pertaining to a given element within a mesh is not solved in isolation. Instead, the sets of equations belonging to all elements of the mesh are collected into a single system of equations describing the global behavior of the problem of interest. Suitable modifications are made to the global system equations to account for boundary conditions. Then, the global system is solved automatically by any of a number of standard methods of digital computing. The ability to solve problems with extremely complex geometry by composing the simple behavior of a multitude of elements is the main innovation of the FEM.
In application to wave propagation physics, the FEM can be applied in several ways. Physically, waves propagate in time and space, so the FEM may be applied directly to this approach. Typically, the FEM is used to model the spatial variation of the wave field, and a separate time-integration or stepping algorithm is used to advance the wave from one instant to the next. The time-integration algorithms in most common use are categorized as “implicit” or “explicit”, referring to whether the solution for the next instant in time is defined as an explicit function of the previous data, or whether an implicit system of equations, relating old to new data, must be solved. Implicit solution approaches are preferred if time increments between old and new data yield savings which are large enough to offset the high cost of solving an implicitly defined system of equations at each increment. Explicit solution approaches are preferred if short time increments are necessary for physical reasons, such as capturing the transient effects of pulse-type waves. Alternatively, and very commonly in wave analysis, time-domain solutions are not sought at all: the time variable is transformed mathematically (typically using a Fourier or Laplace transform) to an auxiliary variable corresponding to the frequency of a steady-state wave. Full understanding of the wave physics for a given problem of interest is derived from solving the solution at multiple frequencies; each solution is independent of the other, unlike time-domain solutions, which depend on the history up to that instant.
Special challenges arise when attempting to obtain a mathematical solution, either exact or approximate, to any problem involving scattering and/or radiation of waves in an unbounded region (that is, in open or free space or in an environment, such as a large body of water, where the region is perceived as essentially unbounded). For example, the problem may comprise radiation of acoustic, elastic or electromagnetic waves. The mathematical solution to such problems must satisfy a so-called “radiation condition,” often known as the Sommerfeld condition, in order for a unique solution to exist. The condition states that all waves “at infinity” are only traveling outward toward infinity, not inward from infinity. Thus, all the energy in the problem resides in the radiated or scattered waves, which are traveling outward after their interaction with an object; no energy is created at infinity. Note that the radiation condition is a condition that exists “at infinity,” not at a finite distance. Much of the history of computational methods for such problems has been focused on how to obtain approximate, numerical solutions that satisfy the radiation condition to an acceptable accuracy, while not being prohibitively expensive. Truncation of the computational model is required at some distance from the source, because computer resources are always limited. Balancing the size of the model with the requirement to model radiation to infinity is the main challenge in selecting where the truncation is made.
The Finite Difference (FD) Time Domain method can be used for large and small problems as well and is generally considered computationally efficient. However, FEM is preferred for realistic engineering problems because FD typically uses stair-stepped rather than conformal approximation of model geometry and “staggers” the field quantities at different locations compromising the accuracy of a solution. The FEM, because of its versatility in handling objects of virtually any geometric shape and material, is generally the method of choice for modeling the finite part of the problem, i.e., the object and, sometimes, a finite part of the open region surrounding the object. It is a continuing challenge to deliver an accurate model for the remainder of the infinite region, including the radiation condition at infinity. Methods for accurately modeling the infinite region generally fall into five classes, with differing applicability and efficacy for time-domain, frequency-domain, implicit and explicit methods of analysis.
A first class uses boundary integral equation methods (BIEM). Here, an integral equation that satisfies the radiation condition exactly can be applied directly on the outer surface of the object. Its advantage, which seems compelling at first sight, is that the infinite exterior domain is replaced, with no loss in physical approximation, by a relatively small surface, which greatly reduces the computational size of the problem. However, it has a severe disadvantage. The matrices in the resulting BIEM discretization of the integral equation are fully populated, making the computational cost prohibitive for large-scale implicit problems. The BIEM can be accelerated through the use of fast-multipole iterative algorithms, but this adds significant complexity to the software. BIEM approaches are much more successful in the frequency domain than the time domain.
A second class uses exact solutions to the wave equation in open regions, often expressed as infinite series of known functions, e.g., wave functions or multipoles. These are joined to the solutions in the finite region in a manner that approximately establishes continuity along, e.g., a closed boundary surrounding the object. This approach suffers from the same disadvantage as the BIEM, namely, that the resulting discretized equations are fully populated, and hence the computational cost is expensive. Wave function solutions are also more successful in the frequency than time domains.
A third class involves constructing an artificial boundary surrounding the object, then applying a so-called absorbing or non-reflecting boundary condition (ABC) that will make the boundary appear as transparent as possible to all outward traveling waves, i.e., the radiated or scattered fields. The primary advantage is that the resulting discretized equations have sparse matrices. Consequently, this property reduces the computational cost dramatically. However, ABC's, which are applied to a boundary at a finite truncation distance, can only approximate the exact radiation condition at infinity. Therefore, spurious (non-realistic) waves are reflected from the artificial boundary, which then propagate throughout the entire finite region, contaminating the entire solution. Although this effect may be mitigated by moving the artificial boundary farther away from the object, this movement only increases the size of the finite region. Consequently, the computational cost once again increases. However, ABC approaches have been successfully applied in both time and frequency domains, for implicit and explicit time-integration schemes.
A fourth alternative to the aforementioned approaches is the use of so-called “infinite elements.” Infinite elements are in fact finite elements that model a semi-infinite sector of space in the same way that finite elements do. To use infinite elements, one first constructs an artificial boundary surrounding the object. However, instead of using ABCs, modal expansions or the BIEM, one defines a single layer of infinite elements around the entire artificial boundary. Because the “infinite” element algorithm imitates the structure of a finite element, the solution algorithm does not change, and the structure of the system of equations is the same as if the infinite elements were finite elements.
Infinite elements have proven extremely useful in solving problems in acoustics. Such uses for infinite elements are described, for example, in the article by D. S. Burnett, “A Three-Dimensional Acoustic Infinite Element Based on a Prolate Spheroidal Multipole Expansion,” J. Acoust. Soc. Am. 96 (1994) 2798-2816 (BURNETT 1994). Acoustic infinite elements are advantageous in that they exhibit both high accuracy and high speed of computation, in frequency- and time-domain approaches to wave problems.
The idea of coupled finite element (FE)/infinite element (IE) approximation for exterior wave propagation problems dates back to the pioneering contributions of Bettess, and Bettess and Zienkiewicz. The works of Astley et al., Cremers et al., Givoli and others recognized the spectral character of the approximation and pointed to the necessity of multipole expansions. Burnett advanced the approach from the practical point of view, by introducing a new, symmetric unconjugated formulation using prolate and oblate spheroidal elements. Contrary to the concept of Perfectly Matched Layer (PML; see below), and other techniques based on Absorbing Boundary Conditions (ABCs), the conjugated element of Astley et al., aims at obtaining the solution in the complete unbounded domain.
A fifth alternative to an absorbing boundary “condition” is an absorbing boundary “layer.” An absorbing boundary layer is a layer used in a model wherein an artificial absorbing material is placed adjacent the edges of the grid independent of the boundary condition. A wave entering the absorbing boundary layer is attenuated by the absorption and decays exponentially such that any reflection is exponentially tiny. This is a numerical abstraction of the absorptive layers used in anechoic chambers to simulate radiation to infinity by minimizing reflections from the chamber walls.
One type of absorbing boundary layer is known as a perfectly matched layer (PML). The key property of a PML that distinguishes it from an ordinary absorbing material is that it is designed so that waves incident upon the PML from a non-PML medium do not reflect at the interface between the PML and the non-PML. This property allows a properly designed PML to strongly absorb outgoing waves from the interior of a computational region while minimizing reflection of the waves back into the interior of the model.
PML was originally formulated by Berenger in 1994 for use with Maxwell's equations for electromagnetic waves, and since that time, there have been many reformulations of PML for both Maxwell's equations and for other wave equations. Berenger's original formulation is known as a split-field PML, because it divides electromagnetic fields into two auxiliary fields in the PML region. A later, more popular formulation, because of its simplicity and efficiency, is called uniaxial PML or UPML (Gedney, 1996), in which the PML is described as an artificial absorbing material, wherein the properties of the PML are different in different directions, e.g., crystals that measure differently or behave differently along two or more axes.
Although both Berenger's formulation and UPML were initially derived by manually constructing the conditions under which incident do not reflect from the PML interface from a homogeneous medium, both formulations were later shown to be equivalent to a more general approach: stretched-coordinate PML (Chew and Weedon, 1994; Teixeira and Chew, 1998). In particular, PMLs were shown to correspond to a coordinate transformation in which one or more coordinates are mapped to complex numbers. More technically, this constitutes an analytic continuation of the wave equation into complex coordinates, replacing propagating, oscillating waves by exponentially decaying waves. This viewpoint allows PMLs to be derived for inhomogeneous media such as waveguides, as well as for other coordinate systems and wave equations.
PML absorption is, in practice, imperfect and depends on the incident angle of a wave. As a wave approaches glancing incidence with respect to the interface between the PML and the non-PML, the attenuation of the wave goes to zero. Consequently, waves sufficiently close to glancing incidence will have round-trip reflections through the PML, thereby introducing artifacts into the solution. Further, where the medium through which the wave propagates is inhomogeneous, without a uniform composition or structure, the PML can fail to absorb and attenuate the wave to create an appropriate model to generate an ultimate solution.
PML absorbing layers with finite elements generally support fast performance and can be applied for non-convex geometries, but require truncation of the model in the computational domain. Additionally, the PML do not absorb all waves in an exact manner since approximate, nonphysical absorptive equations are necessarily used.
Despite the evolution of a variety of FEM modeling technologies and various approaches for providing efficient approximations for wave propagation in unbounded domains, there still exists a significant need for more advanced solutions to support progressively more complex modeling while simultaneously maintaining and improving accuracy, with inclusion and deployment on both desktop and massively parallel computing systems. A need exists for such a computer-implemented apparatus and method for efficiently and cost-effectively modeling wave propagation systems including infinite exteriors based on time-domain explicit finite element calculations wherein it is desirable for the solution approach to exhibit certain characteristics and features.
A preferred solution approach would minimize boundary reflection artifacts in finite-element wave equation computations associated with time-domain explicit finite element calculations or frequency-domain finite element calculations.
Moreover, it would not alter the fundamental algorithmic and data structures of either approach, so that coding and parallelization efficiencies are realized to the maximum potential.
Further, a preferred solution would support time-domain explicit finite element acoustic wave equation computations using a large memory and a high-speed processor, or a multiprocessor massively parallel computing system for execution of associated numerical simulation software.